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Abstract 

We study motion of vortices in arrays of Josephson junctions at zero tem- 
perature where it is controlled by quantum tunneling from one plaquette to 
another. The tunneling process is characterized by a finite time and can be 
slow compared to the superconducting gap (so that rA >> 1). The dissipa- 
tion which accompanies this process arises from rare processes when a vortex 
excites a quasiparticle above the gap while tunneling through a single junc- 
tion. We find that the dissipation is significant even in the case rA >> 1, in 
particular it is not exponentially small in this parameter. We use the calcu- 
lated energy dissipation for the single vortex jump to estimate the physical 
resistance of the whole array. 
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I. INTRODUCTION 



In recent years dynamics of Josephson-j unction arrays has attracted a lot of interest fTj-Q. 
The Josephson-j unction arrays (which are artificially fabricated networks of superconducting 
islands weakly coupled by tunnel junctions) became model systems for the study of quantum 
phase transitions, i.e. transitions occuring at T — > 0. 

The simplest physical picture of the phase transition in a two-dimensional short-ranged 
Josephson-j unction array is the following. The temperature is lower than the bulk tran- 
sition temperature of the islands, so that each individual island is superconducting and is 
characterized by a phase of the superconducting order parameter. The absolute value of 
the order parameter, the superconducting gap A, is the largest energy scale in the problem. 
The phase variable is conjugate to the Cooper pair charge on the island. When the phase is 
well defined, the charge fluctuates and the array is superconducting. That happens in the 
limit where the Josephson energy Ej, associated with the Cooper pair tunneling, is much 
greater than the Coulomb energy Ec-, which determines the electrostatic coupling between 
the islands that tends to localize the charge carriers. In terms of vortices that means that in 
the limit Ej ^ Eq the vortices form the Abrikosov lattice. In the opposite limit, Ec ^ Ej, 
the Coulomb blockade pins Cooper pairs to the islands, so at low temperatures the array is 
insulating. Since in this phase the charge is fixed, the phase variable fluctuates and vortices 
form a superfiuid. 

Both phases were observed by preparing samples with different values of Ec and -E'j [0 . 
The insulating phase exhibits high values of resistance at finite temperatures, which grow 
as T — > 0. The opposite behavior indicates the superconducting phase. The transition can 
also be induced in the same sample by varying magnetic field. The field-induced transition 
can be experimentally observed in arrays and also in granular superconducting films 0. 

Conventional theoretical picture of the superconductor to insulator (S-I) transition sug- 
gested by M. P. A. Fisher is based on duality between vortices and charges. In this picture 
the transition point between the two phases is characterized by finite resistance, which is 
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predicted to have universal value, proportional to the quantum resistance Rq = h/Ae^. 

Experimentally reported values of the transition point resistance, however, while 
being of the same order as the predicted universal value, differ by as much as a factor of 5. 
Moreover, recent experiments show that the superconducting and insulating phases are 
separated by the wide metallic region, characterized by non-zero dissipation. In particular 
it was found that at low temperatures (T < Tq =100 mK) and in a noncommensurate 
magnetic field array resistance becomes temperature independent and remains finite down 
to the lowest temperatures accessible (10 mK). 

The metallic behavior of the arrays can not be described by the usual duality picture, 
since it ignores the presence of dissipation. Two issues have to be addressed. In terms of vor- 
tices, a metal corresponds to a normal liquid, rather than the superfiuid which characterizes 
an insulator. Vortices, however, are interacting bosons and at low temperatures tend to form 
the Bose condensate. Therefore the first question is how can the zero temperature normal 
liquid exist. The second question is what is the origin of dissipation at zero temperature. 

In this paper we will focus on the second question. We consider vortex motion at zero 
temperature where it is controlled by quantum tunneling of single vortices from one plaquette 
to another. It turns out that during the tunneling process a vortex can excite a quasiparticle 
state above the gap with the probability which is not exponentially small in the parameter 
rA >> 1, where A is the superconducting gap on the island and r is the tunneling time. 
The relaxation of the excited quasiparticle then provides the dissipation in the system. 

In order to calculate the matrix element for quasiparticle excitation during vortex tun- 
neling we first solve a simpler quantum-mechanical problem. We consider a particle in a 
quasiclassical potential barrier which is also coupled to a single harmonic oscillator. The 
probability to tunnel through the barrier is given (in the simplest approximation) by the 
WKB approach. The initial state of the whole system (the particle and the oscillator) is 
that before particle tunneling the oscillator was in its ground state. After the tunneling the 
oscillator could remain in its ground state or it could be in one of its excited states; the 
latter case corresponds to dissipation because for any non-zero coupling to environment the 
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oscillator will eventually relax to the ground state. Note that this relaxation can not affect 
the tunneling since it has already happened. For such a problem the dissipation is deter- 
mined by conditional probability of the oscillator excitation (given the fact of the tunneling) 
which we calculate below in Section IV. 

The solution of this quantum mechanical problem can be applied to the case of vortex 
tunneling. When vortex moves the phase on the islands changes. The time derivative 
of the phase acts as an effective field acting on the quasiparticles and thus may result in 
quasiparticle excitations. We note that the processes of excitation of different quasiparticle 
modes are independent, so the result for the total dissipation is given by the sum over all 
modes. 

Besides the calculation of the matrix element, we have to make sure that in the process 
of quasiparticle excitation the energy is conserved. In real arrays the quasiparticle gap A 
is larger that both Ec and Ej, so at zero temperature a single vortex does not have enough 
energy to excite a quasiparticle. However in the vortex liquid dissipation does happen. 
Experimental evidence |^ suggests that the vortex lattice melts easily due to frustration 
(since in the incommensurate magnetic field the vortex lattice does not match the underlying 
array). That means that kinetic energy of vortices in the liquid is rather small, so that the 
liquid retains the short range order. This viewpoint is supported by numerical simulations 
of 2D melting, which show that below the melting point the vortices are mostly in 
large, ordered clusters. When vortices move, these clusters move as a whole and the extra 
momentum due to the interaction with quasiparticles in the island is transfered to the whole 
cluster. The energy of the cluster is much larger then the gap in the quasiparticle spectrum, 
therefore the energy conservation is satisfied. 

The rest of the paper is organized as follows. In Section II we describe Josephson 
junction arrays in terms of phase variables. In the following Section we derive the form 
of the interaction between vortices and quasiparticles. In Section IV we solve a quantum- 
mechanical problem of a particle in a barrier-like potential coupled to a harmonic oscillator 
and in the next Section apply the results to the vortex tunneling. In Sections VI we discuss 
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the vortex lattice melting in the presence of frustration. In Section VII we obtain the array 
resistance due to the dissipation found in Section V. The conclusions follow in Section VIII. 

II. JOSEPHSON-JUNCTION ARRAYS 

The microscopic theory of superconductivity in each individual island will be reviewed in 
the next section; here we shall describe an array of small superconducting islands assuming 
that the amplitude of the order parameter in each island is constant and it is entirely 
controlled by a single phase variable, i.e. we ignore its spatial dependence on the length 
scale of the size a of the islands. This is true when the magnetic field does not penetrate the 
bulk of the islands, which is guaranteed by the condition that the flux through one island is 
less than the flux quantum Ho? < 0o- 

The array Hamiltonian consists in general of three parts. Time variations of the phase 
in each island result in voltage differences between the islands, which are electrostatically 
coupled to each other and to the ground plane. That defines the first part of the Hamilto- 
nian, the electrostatic energy as (1/2) Z^C'jjVi^j) where Cij is proportional to a capacitance 

■i-i 

matrix C^j^ = (2e)^Cj^^. In real experiment [0] the main contribution comes from the junc- 
tion capacitance C. Taking into account also the self-capacitance Cq (the capacitance to the 
ground plane) we approximate Cij by a matrix which only non-zero elements are diagonal 
Cii = Co + 4C and those corresponding to nearest neighbors in the array Cij = —C. The 
junction capacitance defines the energy scale Ec = e^/2C, which is usually referred to as 
charging energy. 

The second part of the array Hamiltonian is the Josephson coupling between the neigh- 
boring islands. The coupling defines the other energy scale in the system Ej. The dissipation 
arises from the coupling of the phase variable to some other degrees of freedom in the array. 
We will denote that part of the Hamiltonian as Hint and will derive its form in the next 
Section, where we consider coupling of the time-dependent phase to the quasiparticles in 
the islands. The Hamiltonian therefore is given by 
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^=2^ ^^i'^ih -EjY. cos(0i - 0,) + H^r^t. (1) 

ij <ij> 

At small Coulomb energy the ground state of Eq. (|I|) corresponds to the constant phase. 
Excitations around that ground state are small spin-wave-like fluctuations and topological 
defects, vortices, where the sum of all gauge-invariant phases around a vortex adds up to 27r. 
In the superconducting state, while the Coulomb energy is small, vortices appear in bound 
pairs (with anti- vortices) . As the Coulomb energy increases, pairs unbound, resulting in 
the transition to the insulating state of the array. In the field-tuned transition, Ec is kept 
constant and vortices are created by the magnetic field. In the superconducting state they 
form a lattice, which melts at the transition point. Melting occurs at Ec which is smaller 
then needed to unbound the vortex-anti-vortex pairs. The density of vortices is controlled 
by magnetic field and remains small even in the liquid phase. The flow of this vortex liquid 
results in finite resistivity. 

Vortices forming this liquid are distinguished by two important features. First, they do 
not have a normal core region, which would be the source of dissipation in homogeneous 
superconductors. Second, since the Coulomb energy (which is a measure for the kinetic 
energy of vortices) is small, the vortex motion at low temperatures is due to quantum 
mechanical tunneling through the cosine potential. 

The tunneling rate can be determined by calculating the instanton action corresponding 
to a vortex moving from one site to the neighboring site. The instanton action was deter- 
mined by several authors [p!0|-[T2| each for slightly different models without dissipation, which 
also differ from Eq. (|l]) by another form of capacitance matrix; note, that theoretical cal- 
culation []T0[ involves also approximation of the cosine Josephson interaction by a piecewise 
parabolic potential, i.e. the Villain's approximation. The results are similar, the instanton 
action is Sinst = cn^J Ej/ Ec, where the number a is of the order of unity and depends on 
a particular model. To determine the resistivity in the array we need to take into account 
also the dissipation (described by the Hint term in Eq. (|I])), which is the main subject of 
this paper. 
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For better comparison with experiments we need to determine the constant a in the 
reahstic model with experimental values of coupling constants. We have repeated the direct 
numerical instanton action calculation for the Hamiltonian Eq. without the dissipation 
term Hint- First we find the phase configuration in the array corresponding to one vortex in 
a particular (arbitrary, but known) plaquette. To do that we set the magnetic field through 
the array so that the total flux is exactly equal to one flux quantum and then minimize the 
energy Eq. Tunneling corresponds to changing the phase configuration to the one with 
the vortex in the neighboring plaquette. In terms of phases the vortex tunneling can be 
described as tunneling of each individual phase in the vortex configuration from it's value 
corresponding to the original position of the vortex to the value corresponding to the final 
position of the vortex. We set these two vortex configurations as the boundary conditions 
for time evolution of individual phases in the array and minimize the action, corresponding 
to the Hamiltonian Eq. (|I]), taking the phases to be functions of imaginary time. 

For array sizes 6x6, 8x8 and 10x10 we determined the value of the coefficient a = 0.7 
(for Co = 0). The phases that change the most during tunneling are the ones in the plaquette 
with the vortex. Therefore even with relatively small array sizes the calculation gives the 
answer that does not change with increasing the size. 

The tunneling rate Fq ~ exp{—Sinst) is then given by 

Fo ~ ^ ^EjEcexpi-0.7^Ej/Ec) (2) 
and provides a measure for the vortex mass. 



III. SUPERCONDUCTIVITY IN A SINGLE ISLAND 

In this section we briefly derive the Hint part of the phase Hamiltonian Eq. (|I]) which 
couples phase fluctuations to quasiparticles. In this derivation we follow the standard mi- 
croscopic description of superconductivity based on the BCS Hamiltonian |T3|,IT^ . 



We start with the BCS Hamiltonian with an effective local, attractive interaction 
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Hbcs = - j d'x^l{x)^M^) - II d'x^i{x)ijUx)^p.^{x)Mx). (3) 

A summation over spins is implied. The order parameter A (a;, r) is introduced by means 

of the Hubbard-Stratonovich transformation. The grand canonical partition function 
Zg = Tr^{exp[—P{H — fiN)]} becomes 

Zg = |/ 1^1^, A*] Texp J drH.ffir)^ | , (4) 
where the effective Hamiltonian is given by 

Heffir) = J d'x{£{x) - f)j Mx) (5) 



-A*{x,T)ij^{x)ij^{x) - A{x,r)ri{x)r^{x) + -|A(x,r)|2}. (6) 

9o 



We can compactify our notations by introducing Nambu spinor [15 



V^=|'^'l- (7) 



The effective Hamiltonian becomes 
Heffir) = f d'x 



(i^fa- A)^ + -|A(x,r)p 



(8) 



where K = —- /i is the kinetic operator and the order parameter is described by 

the matrix A = | A|e^^*^^3fi, |A| and are the absolute value and the phase of the order 
parameter. In our approximation |A| is constant and (/)(r) is a function of time. 

The Hamiltonian is now quadratic in fermion fields and we can formally perform the 
trace over the fermion variables in Eq. (^). The partition function becomes the integral over 
the order parameter 

Zg = Jv[A,A*] exp(-5[A]), (9) 

where the action is 
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S[A] = -Tr\nG-^ + / rfr-|A(r)|2. (10) 



Here G is a 2x2 matrix Green's function in the particle-hole space typical for supercon- 
ductivity which inverse is given by 

G-\x,r;x',T') = 1^--^- Kf3 + A^5{x-x')S{T-T'). (11) 

The action Eq. (|T0|) is a standard BCS action written in the form convenient for the following. 
It does not contain any dissipation as yet, therefore there are no non-local terms discussed 
in Ref. [Q. The dissipation appears after an additional assumption about time dependence 
of the phase variable, namely that while on average it changes slowly, this change occurs 
with rare but large enough jumps, due to the lattice structure of the array. Thus in the 
following we shall not assume that 0(t) is a smooth function of time; such assumption would 
eliminate all dissipation sources in this problem. 

The dependence of the Green's function (and therefore the action) on the phase of the 
order parameter can be displayed through the gauge transformation 

g-\x,T;x',T') = e-^'^^^l'^G-\x,T-x\T')^'^^^l'^, (12) 

where is obtained from Eq. ([TT|) by the replacement — > — t^t's- This transfor- 

OT or 2 OT 



mation shows that a constant (f) contributes nothing to the to the action Eq. (|TT]). 

The fermion contribution to the action can be represented by the path integral over 
Grassman variables 

exp(Tr Ine?-^) = J V[^,] exp{-S^), (13) 

where the fermion action 

Si^ = J drdr' J d^xd^x' iJ^G'^x, r; x', rO^/- (14) 



is explicitly given by 



= J drd^x ip"^ 



^. (15) 



At zero temperature we can write the real time action as 

/O \ 0(p 



The corresponding Hamiltonian in momentum space is given by the 2x2 matrix 

-ek-(p |A| 

|A| tk + ip, 



(16) 



(17) 



where ip = 2~g^' BCS theory (f = and the Hamiltonian can be diagonalised by the 

Bogolyubov transformation, which is just a rotation of the fermion variables. When f ^ 
we still can perform the rotation, but the resulting action will no longer be diagonal due to 
the time dependence of ip. The rotation matrix, which diagonalises the Hamiltonian at each 
moment of time is 

1 / Afc + efc + 

|A| 



y2Xk{Xk + efc + v?) 
After the rotation the Hamiltonian becomes diagonal 

-Afe 



-|A| 
Xk + ek + if, 



(18) 



n 



A. 



(19) 



with the eigenvalue Aa 



'efc + V9)2 + |Ap. The action becomes 



S^ = i dt£k ^\k,t) 



-i 1-7-^ 

dt 2Xk{Xk + ek + v) 



ek + ip -|A| 
|A| efc + v5. 



dip 
'dt 



lik,t), (20) 



where ^{k,t) are the variables in the rotated basis (which for = correspond to Cooper 
pairs) . 

The additional term in Eq. (pO|) appeared due to the time dependence of (p. The diagonal 
part is small compared to the eigenvalues A and can be ignored. The non-diagonal part, 
however, describes a new process : a quasiparticle excitations (the Cooper pairs correspond 
to diagonal part of Eq. (^) ). This is the interaction term which is responsible for dissipa- 
tion. Upon integrating out the fields 7(fc,t) it becomes the interaction part of the action, 
corresponding to the Hint part of the array Hamiltonian Eq. (|l|). 
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S,^, = z f dtd'k^^{k,t) \ f ° ^lik,t), (21) 

Here we have integrated the interaction term by parts, in order to express the resuh in terms 

of the phase fluctuations 0. This brings the extra factor 2Afc from the time dependence of 

7(/c, t). Also we neglected ip in the prefactor which forms the coupling constant gk ~ 
lAI 

— because we consider only adiabatically slow motion. 

2(Afc + e,) 

The phase now can be treated as independent variable, describing a "particle" in the 



periodic potential and coupled to the quasiparticles in the island through the action Eq. (pi|). 
Note that the action Eq. ( ^I]) is diagonal in momentum k, so for each k the quasiparticle 
action is that of a two-level system. Since the probability to excite a quasiparticle is small, 
the phase fluctuations excite only one two-level system at a time, so these excitation processes 
are independent and the total probability can be found as a sum of probabilities to excite each 
individual two-level system. Therefore we can consider the quantum mechanical problem of 
a particle coupled to the two- level system and then integrate the results over k. 

Furthermore, dissipation resulting from exciting a two-level system is not different from 
the one of an oscillator because excitations of the latter to higher levels can be neglected. The 
latter problem has a slightly broader application. Note, however, the important difference 
between this problem and the Caldeira-Legget ||16| model. Here the motion of the particle 
is coupled to a single oscillator with large level spacing, which corresponds to the large 



quasiparticle gap, whereas the Caldeira-Legget jig model is a system coupled to a large 
number of small oscillators, so it's not difficult to excite each one individually and interesting 
physics arises from exciting a large number of them simultaneously. 

Thus we reduced our problem of calculating the probability to excite a quasiparticle 
in an island to a problem of exciting a harmonic oscillator during tunneling. In the next 
section we consider this simpler problem and then in the following section we apply the 
obtained results to the case of coupling to a two-level system and then sum the probability 
over momenta k to obtain the final probability to excite a quasiparticle. 
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IV. SIMPLE MODEL - PARTICLE COUPLED TO HARMONIC OSCILLATOR 



In this section we consider the quantum mechanical problem of a particle coupled to a 
harmonic oscillator and tunneling through some barrier. The Hamiltonian is 

^ = 1^ + ^^''^ + + ^^^^oQ' + 9PQ. (22) 

where p and P are momentum operators of the particle and the oscillator, m and M are their 
respective masses, V{x) is the potential which we assume has a form of a barrier and g is 
the coupling constant. Here we chose the coupling coupling (using the momentum operator 
p rather then the position operator) which has the same form as the one in the action Eq. 
(pT|). We need to obtain the probability to find the oscillator in its first excited state after 
the particle have tunneled through the barrier if before the tunneling the oscillator was in 
the ground state. The coupling g is taken to be small enough so that the oscillator states 
are unchanged. 

We look for the wave function of the system in the form 

^ = Y.^>^{x)\n), (23) 

n 

where \n) denotes oscillator wave functions corresponding to n-th energy level. 
Neglecting the coupling completely we should have 

$0 = vI/o(x) |0), (24) 

which means that the oscillator is in the ground state. The particle wave function under the 
barrier is given in the WKB approximation by 

^o{x) ^ uq{x) exp I - y ^2m{V - E)dx . (25) 

\ Xa / 

In the first order in g we have for the wave function $i = '^o[x)\0) + \E'i(x)|l). The 
Schrodinger equation for the correction \E'i(x) is 



1 

+ V{x) -E + uq 



2m dx"^ 



^i(x) -2(7gio^^o(a;) = 0, (26) 
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where Qio = {1\Q\0) = l/y2MujQ is the oscillator matrix element. The oscillator ground 
state energy ujq/2 is included in the definition of E. It is convenient to express the solution 
in the form \E'i(a;) = ui{x)'$o{x) with the boundary condition ui{xa) = 0, noting that the 
oscillator was in the ground state prior to tunneling. The equation becomes 

^ 7^'/ + u[ - mujo^ui + imgQio = 0. (27) 



2^^ ^ ' 

Compare now first and second derivative terms. The typical particle energy is 

/ V *o 
Q = \ — — , where L ~ (xb — Xa)- The typical time r ~ The ratio — can be es- 

V mL^ Wq 

timated using Eq. (|25|) as , Therefore the second derivate term can be estimated as 

VmV 

^0 1 

— -u" ~ Ml . The first derivative term is simply u[ ~ Mi/L, so that their ratio is 

2\l'o L'^VmV 

\l>o 1 

— -u'l/u[ ~ — j= <^ 1. Therefore we can drop the second derivative term in Eq. (p7|). 
2^Q LymV 

Solving the remaining first order equation we get for the function mi right after the tunneling 



ui{xb) = -imgQio J ^^P I ~ / 



mujo 

=as 



dx. (28) 



2m{V - E) 

Consider now two limiting cases. When uoq <^ (fast transition), the integrand in 
the exponential in Eq. ( P5| ) is small, therefore Ui{xb) ~ —irngQioL so there is no addi- 
tional suppression other than the smallness of g. In the opposite limit ujq ^ Q we can 
evaluate the integral in the exponential noting that the main contribution comes from the 
region near the ending point Xf,- We can then expand the integrand around Xb to get 

. Evaluating the integral is now 



muJo 2mujQ dV 

— ~ . yxb — X, where V = — — 

2m{V - E) V2mV' dx 

straightforward and we get 



x=x, 



v / ^ \ ^ 

ui{xb) = -imgQioL imgQioL — , (29) 

which means that in the limit of slow tunneling Mi (a;;,) is indeed small in Q/ujq but only as 
a power law. 

The probability of exciting the oscillator is proportional to the square of ui{xb) and is 
given by 
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V^gVL'Q\J^)\ (30) 

The average energy W dissipated in one jump is equal to the energy needed to excite the 
oscillator [uoq) times the transition probability Eq. (pO|) 

W ^ {gmLQ^^f (31) 



V. TUNNELING OF PARTICLE COUPLED TO QUASIPARTICLE SYSTEM 

We now return to the full problem, formulated in Section ||, namely to the action Eq. 
(0). As we have already mentioned, for each value of k we can treat the quasiparticle 
system as a two-level system and then sum over all possible k. Again we treat the phase 
variable on the island as a coordinate of a "particle", which tunnels through some barrier. 
The result for this case is exactly the same as for the case of the oscillator since we have 
neglected excitations of higher levels. In the case of a two-level system there are no higher 
levels at all and the result Eq. ( pT]) is the full answer, in which we have to substitute the 
corresponding matrix element for Qiq and the value of the gap for ujq. 

Therefore in order to apply our results to the case of the phase coupled to the quasiparticle 
system Eq. (ETI) we only need to rewrite it in a Hamiltonian form equivalent to Eq. 



In the Hamiltonian formalism the derivative — in Eq. (^) is replaced by the momentum 
operator p/m [m is the "particle" mass) leading to the effective interaction constant 

9k = — T^TT , 7- (32) 

The energy Q is now the typical frequency of the phase variation and is equal to the inverse 



tunneling time, ~ y^EjEc- 

The probability to excite the quasiparticle at each k then follows from Eq. (^Ol). The 
oscillator frequency Uq is now exchanged for the quasiparticle gap, which is 2Afe at the same 
k. The matrix element corresponding to Qio is just 1. The probability then is 
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(33) 



Here we used the factor of 27r for the effective length L, which is by how much the phase can 
be changed. The dissipation contribution for each k follows by multiplying the probability 
be the energy gap 



where A/" is the number of particles on the island. This defines the energy transfered to the 
quasiparticle system during tunneling and therefore the dissipation in a single vortex jump 
between two neighboring plaquettes. 



At low magnetic fields vortices form a lattice that melts at higher fields. Because melt- 
ing is due to the competition between kinetic and interaction energies, it happens when 
the two are parametrically equal. However in 2D the liquid retains short-range order and 
the interaction energy loss in melting is numerically much smaller as described by a small 
Lindemann number. In arrays the vortex lattice is frustrated by the incommensurability 
with the underlying array structure. This effect reduces the ratio of kinetic energy over 
interaction energy at melting even further. 

In the absence of a microscopic theory of melting we use the phenomenological Linde- 
mann criterion, which describes the melting in terms of elastic constants of the vortex lattice. 
In the array system these constants are renormalized by frustration. To estimate this effect 
we analyze the experiment on thermal melting [^], in which the effect of incommensurability 
on the transition temperature was studied in detail. We emphasize that these measurements 
were performed on array systems, which are different from the ones discussed throughout 




(34) 



Integrating this expression over momentum we get 




(35) 



VI. VORTEX LIQUID 
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this paper. Here we use these experiments to obtain estimates of the renormahzation of the 
elastic constants of the vortex lattice and then use this renormahzation to describe quantum 
melting. Our observations however are general and therefore are applicable to the quantum 
systems of interest J^. 

We need to estimate how frustration renormalizes the elastic constants. Therefore we 
estimate the interaction energy in the experimental system [H, which translates into un- 
renormalised values of elastic constants, leading to an estimate of unrenormalised melting 
temperature Tmo- Comparing T^o with the experimentally observed Tm we find the frustra- 
tion factor. 

To estimate the interaction energy we relate the superfluid density ps to the observed 
magnetization M. Energy and current in the Josephson junction in magnetic field can be 
written as 

rn ^ T T T • / 2eAa\ 

E = — Ji COSV9 , J = Jism [ if -\ — - — (36) 
2e \ he J 

where A is the vector potential and a is the distance between islands. We relate the Josephson 

energy Ej = ^Ji to the magnetization M at low fields, where response is linear. We use 

the relation j = psA between the supercurrent and the superfluid density ps, and express M 

via the current. Assuming for simplicity circular geometry the magnetization is given by 

L 

M = — J 27irdr{j x f) (37) 



Taking the integral we relate the magnetization to the magnetic field H and the superfluid 
density ps 

M = YcP^\l'H (38) 

We use Eq. (^) to deduce the value of ps from the data |P; for its zero temperature value 
we get Ps(0) = 6.49 x 10^^. The superfluid density is a function of reduced temperature, 
PsiT) = ps(0)r. Since we need the superfluid density at the true melting temperature Tm, 
T is deflned from the shift in melting temperature due to frustration and was measured to 
be r = 0.01. 
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Comparing the supercurrent equation with Eq. (|36|) we obtain the Josephson energy (and 
the magnitude of current) 

E = ^/r = gp.(T.). (39) 

We can now estimate the interaction energy to be ~ 1.7 x 10^ K. Because the melting 
temperature is of the order of 1 K, its ratio to the interaction energy Tc/E = C is estimated 
as C = 10"^ 

We now compare Lindemann criterion for thermal and quantum melting. For the thermal 
melting considered above we have 

J2 



ipp) ~ T / ^ = ai. (40) 



where the integral is over the Brillouin zone. Here we have also substituted ceei?^ for the 
actual dispersion law. This rough estimate will be sufficient for our purposes. Assuming 
that frustration renormalizes the elastic constant by cge ~ ups, we get the renormalisation 
factor K, ~ T / psa\ C^ja^l. Taking for the Lindemann parameter the usual value ~ 0.1 
we get K = 0.01. 

For quantum melting we conjecture that because the renormalisation k is due to frus- 
tration (induced by the array) the elastic constant will be renormalized by the same factor. 



1^ + CmT V C66 

Using the renormalisation factor k we get Ecj Ps ~ kxi\. Thus we expect that the quantum 
melting happens at kinetic energies which are at least three orders of magnitude smaller 
than the interaction energy. 

VII. DISSIPATION EFFECTS ON VORTEX MOTION 

Before we turn to the estimate of array resistance, we have to address the question of en- 
ergy conservation. The quasiclassical calculation of the probability to excite a quasiparticle 
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considered above assumed implicitly that the energy of the vortex is large enough. Quanti- 
tatively that means that the vortex energy should at least be larger than the quasiparticle 
gap, otherwise the vortex would lack the energy to excite the quasiparticle. 

The energy acquired by a vortex driven by the Lorentz force is proportional to the applied 
current so at very low currents it would not be sufficient to excite quasiparticles above the 
finite gap. For larger currents the probability to excite a quasiparticles is constant and is 
given by Eq. ( ^3l) so at these currents vortex dissipation is linear in its velocity leading 
to Ohmic conductivity. Here we estimate smallest currents, jo, at which the conductivity 
remains Ohmic. 

There are two effects that make jo very small. First, a vortex makes many jumps between 
consecutive emission processes and accumulates energy. Second, vortex liquid is incompress- 
ible and retains short range order in a broad range of fields above the melting point, so each 
emission process slows down not an individual vortex but a large number of them. This 
effect is similar to the Mossbauer effect in crystals, but here the momentum can not be 
transfered to the whole number of vortices since there is no long-range order. The effect is 
difficult to describe quantitatively, due to the absence of a theory of a strongly correlated 
vortex hquid at T = 0. We shall attempt therefore only to show that the effect is indeed 
large. 

To see this effect and to estimate the correlation length in the liquid state we perform 
the following calculation, similar to the calculation of the Debay- Waller factor. The idea is 
that when the momentum is being transfered to the liquid as a whole the quantum state of 
the liquid does not change, so it is described by the same wave function after the transfer 
as before. The amplitude of such process is 



where N is the number of vortices, p is the transferred momentum and \l'(xi...XAr) is the 
macroscopic wave function of the vortex liquid. The factor 6{xi) singles out the coordinate of 
the island,where the quasiparticle was excited. When the vortex interacts with the particular 




(42) 
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island its coordinate becomes fixed and therefore we do not need to integrate over it. 

Tfie amplitude A is a function of momentum p, vortex number and interaction strength. 
In a true liquid, where no order is present, the wave function "^Ixi-.-Xn) depends only on 
the relative coordinates of vortices and therefore A = 0. When vortices are organized in 
clusters at short distances, the correlations decay exponentially like exp(— x/^), where the 
correlation length ^ defines cluster size. Therefore for a system of the finite size L the 
amplitude A is of the order of unity when L ~ ^, but when the system becomes large, so 
that L ^ ^, then the amplitude is exponentially small A ~ exp(— L/^). 

Because the vortex-vortex interaction is proportional to the logarithm of distance be- 
tween vortices the vortex liquid can be approximated by the two-dimensional Coulomb gas. 
The exact wave-functions of the Coulomb gas are unknown. However there exists a three- 
body Hamiltonian with interaction that is Coulomb for long distances, while different (and 
three-body) for short distances. The ground-state wave function for this Hamiltonian is 
known [|17|, and it was argued that one can use this known wave function to estimate the 
properties of the Coulomb gas. The amplitude A then can be estimated numerically. Due 
to the limited computer availability we performed the calculation for relatively small arrays 
of up to 20 particles in the circular geometry and small interaction parameters a = 0.5^2.5 
(for comparison, the melting point is at a = 30 |^). But even being that far from the 
melting point we could see the momentum dependence of the amplitude A as described 
above. Our calculation allows to estimate the correlation length ^ > 5 in the units of lattice 
spacing, which would correspond to the cluster size of up to 20-30 vortices and increasing 
as we increase a. 

We now estimate the average energy Ed of a moving cluster due to the external current. 
When the current is small the force acting on a vortex is given hj F = J^q/qc, where a in 
the lattice spacing, J is the current per junction and $o is the flux quantum. The energy 
acquired by the vortex after one tunneling jump to the neighboring plaquette is equal to 
the force times the lattice spacing E = Fa = J^q/c. Since the probability to excite a 
quasiparticle is small the excitation is a rare process and the average number of jumps the 
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vortex makes before it excites a quasiparticle is inverse probability. Therefore at the moment 
of quasiparticle excitation the vortex would have the energy E = J^q/cV. Multiplying this 
energy by the number of particles in the cluster we obtain an estimate of the average energy 
of the cluster at the moment of quasiparticle excitation Ed = N^J^o/cP. 

The probability can be estimated from Eq. (^) using the experimental |^ values for 
Ec and |A|. We get V ~ 0.001. Therefore the cluster energy can be estimated as 
Eel = N^J ■ 1 X 10^ K, where the current is measured in nanoampers. 

The system has linear response when the average cluster energy is larger than the gap 
|A| ^ 2 K. For very small currents, when Ec < |A| the system would exhibit non-linear 
current-voltage curves. For a cluster size A'^ = 100 the current value where this non-linearity 
would be observable can be estimated by setting E^ equal to the gap and using the above 
estimate for the cluster energy. We get Jq = 0.1 x 10^^ nA. The currents used in experiment 
are of the order of J ~ 0.1 nA. Therefore it is likely that in the experimentally observable 
case the array is in the (pseudo) linear regime. 

Thus the coupling to quasiparticles results in dissipation described by linear response in 
contrast to the dissipation due to coupling to the acoustic phase modes (spin waves) |]TB|,|TU[. 
The effective action obtained in Ref. [|r^ has a dissipation term which is proportional to 
tu^ Inuj. For the slow vortex motion we consider (due to very small currents discussed above) 
this term is small compared to the linear term Eq. (35) which is determined by the energy 
scale Ej rather than the frequency. For larger currents the situation is different and the two 
effects become comparable. 

For the (pseudo) linear response regime we now consider how the vortex motion is affected 
by the dissipation. Our goal is to obtain an expression for the total array resistance which 
arises due to dissipation Eq. (p5D, therefore we now consider a macroscopic equation of 
vortex motion, averaged over the whole array following Ref. |^0|. Note that after averaging 



this equation (Eq. (|43|)) does not describe the microscopic couphng between vortices and 
quasiparticles, which in general is non-linear. 

Under the influence of the driving current J a vortex is moving in a direction perpendic- 
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ular to the current flow: 

V^^x + WV^^x = ^, (43) 

where x is the vortex position $q is the flux quantum, Fq is the tunnehng rate, deflned in Eq. 
(0), which provides the measure for the vortex mass and a is the lattice spacing. The lattice 
potential was taken into account when we calculated the dissipation and the tunneling rate, 
therefore it does not appear in Eq. (|43|) . 

If the driving current is constant then the vortices move with constant velocity. If in 
the Cartesian coordinate system the current flows along y axis, then vortices move along x 
(magnetic field is along z axis, perpendicular to the xy plane of the array) and their velocity 
is 

V, = ^r„. (44) 

The potential difference caused by the time dependent phase is 

t/ = -#. (45) 

2edt ^ ^ 

When one vortex is moving across the array (in time t = d/vx, where d is the size of the 
array) the phase changes by 2tt. To obtain the total voltage the effect of one vortex should 
be multiplied by their number = Bd'^/ 

f/ = vr^ ^ J (46) 

The current per junction J can be obtained from the total current J as J = la/d (assuming 
a square array). The coefficient of proportionality between the voltage U and the current / 
is the array resistance 



21 



VIII. CONCLUSIONS 



We have considered the vortex tunnehng in Josephson-j unction arrays at zero temper- 
ature. Using the simple quantum-mechanical analogy we have showed that such tunneling 
is accompanied by small dissipation due to quasiparticle excitations in the superconducting 
islands. Even in the presence of the large quasiparticle gap the probability of such excita- 
tions is found to contain only power-law smallness in the (small) ratio of the characteristic 
tunneling frequency to the gap. 

Our main result is the resistance of the vortex liquid Eq. (^7]), which is due to this 
dissipation. The result is valid for not too small driving currents (for experimental setup 
we estimate J > Jq = 0.1 x 10^'^ i^A), for which the system is in the (pseudo) linear 
response regime. Our results provide the quantum-mechanical mechanism of dissipation in 
Josephson-j unction arrays. The resistance Eq. (^T]) is independent of temperature (due to 
it's quantum-mechanical origin) in agreement with the experiment [0. 

Our argument is applicable to the vortex liquid just above the melting point, where the 
vortex liquid retains short-range crystalline order. Then the vortices are strongly interacting 
and the external momentum can be transfered to a large number of vortices. This is similar 
to the Mossbauer effect in crystals where the external momentum is transferred to the whole 
crystal. In the vortex liquid such a transfer is impossible due to the absence of long-range 
order, but the transfer to a finite size cluster remains possible. Such a cluster involves a 
large number of vortices, therefore its energy is much more than the gap in the quasiparticle 
spectrum, allowing excitations above the gap. Estimating the minimal cluster size from our 
numerical data we got that in experimental conditions the energy contained in such cluster 
is always larger than the gap; it would get smaller than the gap only for very small currents 
J < Jo. 

However these estimates depend crucially on the structure of the strongly correlated 
vortex liquid formed when the vortex lattice melts. We have argued that the problem 
is exacerbated by the frustration imposed on the vortex lattice by the underlying array 
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structure. The frustration reduces even further the kinetic energy needed for the melting. 
The self-consistent theoretical description of the normal liquid of vortices, however, remains 
to be an unresolved question. In particular it is not clear whether the existence of the normal 
liquid is due to the dissipation effects or it is in fact possible to form a normal liquid in the 
absence of dissipation. 

One of the possible descriptions of the strongly correlated vortex liquid is a dilute gas 
of dislocations in the vortex crystal similar to the hexatic phase appearing in 2D thermal 



melting . In this phase the vortex flow is due to the motion of dislocations. The motion of 
each single dislocation transfers the whole row of vortices across the system. Here the number 
of moving vortices scales with the system size so in the thermodynamic limit the combined 
energy of these vortices becomes infinite and linear response persists to zero currents. 

The detailed description of the strongly correlated vortex liquid is a subject of future 
work. 
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